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ABSTRACT 

We have applied an eccentric accretion disc theory in simphfied form to the case 
of an accretion disc in a binary system, where the disc contains the 3:1 Lindblad 
resonance. This is relevant to the case of superhumps in SU Ursae Majoris cataclysmic 
variables and other systems, where it is thought that this resonance leads to growth of 
eccentricity and a modulation in the light curve due to the interaction of a precessing 
eccentric disc with tidal stresses. A single differential equation is formulated which 
describes the propagation, resonant excitation and viscous damping of eccentricity. 
The theory is first worked out in the simple case of a narrow ring and leads to the 
conclusion that the eccentricity distribution is locally suppressed by the presence of the 
resonance, creating a dip in the eccentricity at the resonant radius. Application of this 
theory to the superhump case confirms this conclusion and produces a more accurate 
expression for the precession rate of the disc than has been previously accomplished 
with simple dynamical estimates. 
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A superhump is a periodic variation in the luminosity of an 
accreting binary system, with a period close to the orbital 
period of the binary. Superhumps were first discovered in 
the SU Ursae Majoris class of dwarf novae, where they oc- 
cur during superoutbursts which are longer and more pro- 
nounced than common outbursts. These usually show so- 
called 'normal' or positive superhumps, where the super- 
hump period is slightly longer than the orbital period. Neg- 
ative superhumps, where the superhump period is shorter 
than the orbital period, have also been observed. Some 
other binary systems, such as V603 Aquilae (Patterson & 
Richman 1991), have been described as permanent super- 
hump systems, where the superhumps do not occur in out- 
bursts but are instead variations in the normal lightcurve of 
the systems. Superhumps have also been observed in low- 
mass X-ray binaries (LMXBs), for example KV UMa/XTE 
Jl 118+480 (Uemura et al 2002). 

The most probable explanation for superhumps in the 
SU UMa systems is in terms of the effect of tidal stresses on 
a precessing, eccentric accretion disc, a model first proposed 
for the SU UMa systems by Whitehurst (1988) and investi- 
gated by Hirose & Osaki (1990) and Lubow (1991). This last 
paper derived an analytic expression for an instability at the 
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3:1 Lindblad resonance that was first noted by Whitehurst, 
which leads to eccentricity growth in a previously circular 
disc. Haswell et al. (2001) proposed a similar explanation 
for the superhumps in LMXBs, in which the precessing disc 
leads to a modulation in the light from reprocessed X-rays 
from the central source. 

Observational support for the idea of an eccentric disc 
first came from eclipse maps of dwarf novae (Hessman et 
al. 1992; Rolfe, Haswell & Patterson 2000) and has recently 
been added to by the spectroscopic study of Zhao et al. 
(2005). Many numerical models of eccentric discs in cata- 
clysmic variables have also been studied. Smoothed-particle 
hydrodynamics codes (e.g. Hirose & Osaki 1990; Lubow 
1992; Murray 1998) in general produce eccentric discs and 
superhumps, although some Eulerian codes do not (e.g. 
Heemskerk 1994; Kornet & Rozyczka 2000). 

Theoretically, the precession rate of an eccentric disc 
has traditionally been estimated using the dynamical pre- 
cession frequency ojdyn of a particle due to tidal forces. This 
is given by the expression 



t^dyn 



(l+q)i/2 



1/2 



"3/2 



(1) 



as a fraction of the orbital frequency of the binary Clh, where 
a is the semi-major axis of the binary and fcgyj is a Laplace 
coefficient (see equation (1251 below). Observationally, the 
precession rate is usually expressed in terms of the super- 
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hump period excess e, defined as 



gravitational potential <1>, the equations of motion are 



fib 



(2) 



Mass ratio and precession rate should thus be roughly pro- 
portional, since equation (|5J can be approximated as 



+ 



(3) 



for the low values of q found in SU UMa stars (Patterson et 
al (2005) found q < 0.35). 

Observations indeed show an approximately linear re- 
lation between q and e, and for example Patterson (2001) 
found 



e = (0.216 ±0.018)(7. 



(4) 



As pointed out by Murray (2000), although equation Q 
reproduces the correct behaviour of disc precession period 
with mass ratio, it does not accurately reproduce the con- 
stant of proportionality, which is about 0.4 if the dynamical 
precession is evaluated at the 3:1 Lindblad resonance itself. 
To obtain a more accurate expression, we therefore need a 
fluid dynamical treatment of the disc. The overestimation 
may be due to a subtractive term in the relation due to 
pressure effects, as argued by Murray (2000) using the anal- 
ysis of Lubow (1992). Alternatively, Patterson (2001) noted 
that the correct proportionality can be reproduced simply 
by changing the radius at which equation Q is evaluated 
(with (r) ~ 0.37a to match the observations) and suggested 
that the precession rate may be determined by the effect 
of the tidal potential on a larger disc annulus rather than 
just the resonance. The actual precession would then be an 
average in some way over the disc. 

The aim of this paper is to construct a more complete 
description of the behaviour of eccentricity in the accretion 
discs of SU UMa systems, using the theory of eccentric discs 
derived by Ogilvie (2001). The following section details the 
derivation of an equation for the eccentricity in a 2D disc, 
including the effects of the tidal precession due to the binary 
companion, eccentricity excitation at the 3:1 Lindblad reso- 
nance and damping of eccentricity due to viscosity. Section 
121 then describes solutions of this equation in the form of 
normal modes for both a narrow ring and a full disc, which 
give values for the disc precession rate and the growth rate 
of eccentricity. These results are summarised and discussed 
in section 21 
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These equations are solved by treating each quantity as 
a function of a parameter e, which is a characteristic value 
of the disc semi-thickness H/r and is assumed to be small. 
Quantities can be expanded in terms of e as power series 



X = xo + X2e + ■ 



(9) 



and the terms at each order in the equations equated. 

In the basic state of the disc, we assume that the disc 
is steady, axisymmetric, non-self-gravitating and has no ra- 
dial motion. We also assume vertical hydrostatic equilib- 
rium (valid in some sense if the 2D quantities are treated as 
3D quantities integrated vertically through the disc), which 
leads to a scaling p{r) ~ ^ p{r) between the pressure and 
density and therefore a factor of in the expansion for the 
pressure. Simplifying equations for this basic state, 

we obtain 
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for the angular velocity fl, where $0 is the potential for 
a point mass and $2 is an axisymmetric correction due to 
the binary companion, taken to be of O(e^) since it is small 
compared to the primary potential. 

We then perturb the fluid equations, introducing a lin- 
ear, adiabatic perturbation x'{r)e~^'^ for each quantity. The 
perturbations satisfy the equations 
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We expand the perturbed quantities in the same way as the 
basic state: 



2 THE ECCENTRICITY EQUATION FOR A 
2D DISC 

2.1 Derivation of the eccentricity equation 

To derive an equation for the eccentricity in a disc, we first 
consider the unperturbed disc, using a simple two dimen- 
sional, inviscid model. In terms of the polar velocity compo- 
nents {u,v), pressure p, density p, adiabatic exponent 7 and 



'1 ' 2 , 

Xo+X2£ +■ 



(16) 



Since we are interested the evolution of the disc on time- 
scales slower than the orbital time-scale, we also introduce 
a slow time variable r, related to the time coordinate as 
r = €^t. We do not need to expand the potential, since, for 
a circular binary system over period of time, the average 
potential is axisymmetric. 

At lowest order, equations 11211 and I13II yield the rela- 
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tions 

—iQoUQ — 2^1ov'q — 
—iiiov'o + ^D,ou'q = 0. 

Therefore 2v'o — —iuq, and we can set 
u'o = irD,oE(r,T) 

By comparison with the expressions for velocities in an el- 
liptical Keplerian orbit, we see that, for small eccentricities, 
this describes a Keplerian ellipse, and the function _E(r, r) 
is the complex eccentricity E = ee'" of the orbit. This solu- 
tion arises because pressure and tidal effects do not enter the 
equation at this order, and we obtain the expected solution 
for a free particle in a Keplerian potential. 

By taking an appropriate linear combination of equa- 
tions il'2l and llljl l at O(e^), y.2 and v'2 can be eliminated 
and an equation for the eccentricity derived. We use equa- 
tions 1141 1 and 1151 1 at lowest order to eliminate p'o and p'o, 
and equation IIH to ehminate Q2 in favour of the tidal po- 
tential. Putting this together and removing the scaling fac- 
tors, we obtain a partial differential equation for the complex 
eccentricity 
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The form of this equation without the tidal potential al- 
ready sheds some light on the nature of the various solutions. 
We can Fourier analyse in time by setting E{r, t) = i5(r-)e"^'. 
The eigenvalue ui is determined by the Sturm-Liouville equa- 
tion 
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Solutions of this equation with appropriate boundary con- 
ditions lead to a series of normal modes in the disc. Each 
mode is a radial eccentricity distribution which precesses at 
a frequency given by its eigenvalue. 



2.2 The tidal potential 

To make further progress, we require the tidal potential due 
to a binary companion. As we are interested in long-term 
evolution, we can use the Gauss averaging method (e.g Mur- 
ray & Dermott 1999) where the companion is treated as a 
smeared-out ring of material along its orbit, with the den- 
sity proportional to the time taken to traverse the element. 
For SU UMa systems the orbit of the companion is circular, 
and thus the potential is given by the expression 
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(23) 



for a companion of mass M2 at polar coordinates (a,(;/>b). 
This integral can be written in a shortened form, using a 
Laplace coefficient from celestial mechanics, as 



(24) 



where the Laplace coefficient fe^^^ is defined as 
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This is then inserted into equation 12H . After some manip- 
ulation and simplification using identities for the derivatives 
of Laplace coefficients, it becomes 
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where q = M2/M1 is the mass ratio of the binary. The 
new term is also proportional to E and simply modifies the 
Sturm-Liouville eigenvalue problem of the previous section, 
leading to a new eigenvalue and corresponding eigenfunc- 
tion. 



2.3 Generation and damping of eccentricity 

The Sturm-Liouville problem of section |5| describes how an 
eccentric disc evolves, but does not consider how the eccen- 
tricity in the disc is produced or describe whether it can 
be damped. Eccentricity in the discs of SU UMa systems 
is thought to be due to an instability at the 3:1 Lindblad 
resonance, and an analytical expression for the eccentricity 
growth rate ^ was found by Lubow (1991): 



C = 2.08Cnb<7 



(27) 



where the binary has an angular velocity Q,h and C is a 
correction factor dependent on the size of the ring. We in- 
sert this growth term into the equation describing the ec- 
centricity as a delta- function at the resonance point, so the 
expression describes the rate at which eccentricity is created 
by the resonance and the fluid dynamics of the eccentricity 
equation itself describe how it propagates through the disc. 
The correction factor C, which describes the effect of finite 
disc width on the growth rate, is therefore 1 and the de- 
nominator is replaced by a delta-function in the numerator. 
Thus the modified form of equation 12611 is 
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(28) 
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When analysed as an eigenvalue problem, the eigenvalue 
is now complex, with a real part describing the precession 
frequency of the disc, and an imaginary part which deter- 
mines whether the mode grows or decays in time. To analyse 
the eccentricity growth more completely, we introduce a vis- 
cosity to the disc, which would be expected to compete with 
the resonant growth and cause damping of the eccentricity. 
Ogilvie (2001) demonstrated that a Navier-Stokes shear vis- 
cosity alone results in an overstability of the 2D disc, and 
since we wish to consider damping of eccentricity, we work 
with a bulk viscosity rather than a shear term. Thus we 
introduce a bulk viscosity described by a Shakura-Sunyaev 
(1973) a-parameter. The equations of motion of the disc 
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JSJ-® become 

du du V du 
dt dr r d(j) r 
dv dv V dv uv 1 dp 1 9<3? 1 dT ,„„, 
at or r o4> r rp ocp r o<p rp ocp 
where the quantity T describing the bulk viscosity is given 

by 

Since this viscous term is proportional to the pressure, we 
scale it to the same order 0{£^) and the lowest order equa- 
tions are unaffected. Therefore the perturbed flow can still 
be described as an eccentricity, now governed by a new ec- 
centricity equation. Repeating the analysis that led to equa- 
tion II21II . we find that the bulk viscosity has the simple effect 
of changing the adiabatic index 7 to a complex number; 



7^7 — IQb. 



(32) 



Therefore the full equation describing generation, damping 
and dynamics of eccentricity in the disc is 
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Splitting this into its constituent terms, we can interpret 
the effect of each on the dynamics of eccentricity. The term 
proportional to 7 has the form of a dispersive wave equation 
like the Schrodinger equation, allowing the eccentricity to 
propagate through the disc as a dispersive wave. By analogy 
with the integrated probability density in the Schrodinger 
equation, this dispersion has a conserved quantity 

d 



dt 



-rpn\E\^dr = 



(34) 



with the boundary condition assumed below {pdE/dr = 
at the boundaries). The dynamical term simply describes a 
precession of the eccentricity which also satisfies equation 
13411 . The last term describes the eccentricity growth, and 
the term proportional to the bulk viscosity describes diffu- 
sion of eccentricity through the disc. These terms now break 
the conservation relation of equation 13411 . and we have the 
more general relation 
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The quantity on the left hand side of equations 1341 - 13511 is 
the angular momentum deficit of the disc, which expresses 
the difference between the angular momentum of elliptical 
and circular orbits with the same energy. 

As before, we can Fourier analyse in time to obtain 
an eigenvalue problem similar to a Sturm-Liouville prob- 
lem, but with complex eigenvalues. This still has a series 
of modes, which can be characterised by their number of 
'nodes' or points where the phase of the eccentricity passes 
sharply through tt. A useful expression for this eigenvalue 
can be derived by multiplying equation 1331 by r^pE* and 



integrating over the radial extent of the disc. The term in the 
eccentricity gradient can be integrated by parts and, with 
the boundary conditions we use below, the boundary term 
vanishes and we obtain expressions containing integrals of 
real quantities only. Therefore we find the following expres- 
sions for the real and imaginary parts of the eigenvalue: 



»(lj) = 
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(36) 
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These expressions show that there is a separation of ef- 
fects into those which affect the precession rate, described 
by the real part of the eigenvalue, and those which affect 
the growth rate, which is determined by the imaginary part. 
The real part has three separate contributions. The first 
two are due to pressure, and if the pressure gradient is neg- 
ative (as in the disc models we used) both of these terms 
are negative and lead to a retrograde precession. The poten- 
tial term is positive and induces a prograde precession, so 
these three terms will compete to determine the direction 
of overall disc precession. The two terms in the imaginary 
part behave as expected, with the resonant forcing acting 
to increase eccentricity and the bulk viscosity damping it (a 
negative eigenvalue implies a growing eccentricity). Equa- 
tion 1371 is simply another form of the general relation for 
eccentricity growth and decay expressed by equation 1351 . 
in terms of a growth-rate eigenvalue. 



3 THE SOLUTIONS OF THE ECCENTRICITY 
EQUATIONS 

3.1 Behaviour of eccentricity in a narrow ring 

The overall equation for the eccentricity can be simplified by 
considering a narrow ring centred on the 3:1 Lindblad reso- 
nance. In this case, the pressure and density can be assumed 
to be approximately constant and the tidal term due to the 
companion can be ignored as it causes a constant precession, 
leaving only the eccentricity gradient and the resonant term. 
Equation 1331 therefore greatly simplifies to the form 



dE 
'dt 



1(7 — iab)p d^E 
2^ 'dx^ 



+ ^ES{x) 



(38) 



where x — r — rras- This is helpful because the equation now 
has an analytical solution and the behaviour of the eccen- 
tricity can be examined directly. The parameter ^ can be 
varied to study the effect of changing the resonant strength. 
Equation 1381 can be solved as an eigenvalue problem with 
an eigenfrequency lj as before. We solve the equations in a 
ring extending between x — —a and x — b, with the bound- 
ary conditions that the eccentricity gradient is zero at the 
edges of the ring. The actual numerical value of the eccen- 
tricity is multiplied by an arbitrary constant, and we scale 
by setting eccentricity at the inner edge —a to 1; since an 
absolute eccentricity >1 describes an open orbit, eccentrici- 
ties greater than 1 imply that the eccentricity at the inner 
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edge must be low. The solutions are 



E{x) = cos k{x + a) 
cos ka 



E(x) 



cos kb 



cos k{x — b), < r ^ b, 



(39) 



The wavenumber k is complex and satisfies the equation 



fc(tan(fca) +tan(fc6)) 



iW 



p(7-iab) (a + 6)(7-iQb) 



(40) 



where we have introduced the dimensionless parameter W 
to characterise the strength of the resonance. The eigenfre- 
quency is then given by 



(7-iab)p . 
2pQ 



(41) 



With no resonance, the eccentricity distributions are simple 
cosine waves. When the resonance is added a cusp devel- 
ops at the resonance point in the centre of the ring, with 
a discontinuity in the eccentricity gradient. The behaviour 
of this cusp as the resonant strength increases is dependent 
on the behaviour of equation H4()|l which can be expanded 
and solved analytically for small and large values of W . For 
<C 1, fc or tan(fca) (and therefore both for the lowest 
order solution) must be small and therefore we can write 
tan(fca) ~ ka. Therefore the wavenumber is approximately 



p(7 - iab)(a + b) ' 



(42) 



with a phase depending on the ratio of 7 to a. Equation 
14111 at lowest order in k gives an approximately imaginary 
eigenvalue 



a + b' 



(43) 



in agreement with the growth rates derived for narrow rings 
by Lubow (1991). The following term in the series gives a 
more accurate expression for the growth rate 



9(c^) 
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hb 3p{-/2 + al){a + b)3' 
and an expression for the precession frequency of the disc 
2S,^n-fp{a^ + b^) 



3p{j^ + al)(a + b) 



(45) 



which is purely retrograde and pressure-driven, as the tidal 
potential of the companion has been ignored. 

For a very strong resonance with ^ 1, the dominant 
terms are tan(fca) or tan(fcfe), which can become arbitrarily 
large as ka or kb approaches (n -I- l/2)n. This leads to a 
lowest-order solution in which (assuming a > b) 



TT 

2a 



(46) 



with the eccentricity strongly suppressed at the resonant 
radius, such as shown in Fig. Q for W = 225. As the res- 
onance is made stronger the cusp becomes deeper and the 
wavenumber closer to the value where the eccentricity has a 
zero at the resonant point. The growth rate converges on a 
constant, since the wavenumber is now approaching a single 
value, and is given by 



3(c.) 



8npa2 



(47) 




Figure 1. The eccentricity of narrow rings with different resonant 
strengths, showing the development of the cusp as W increases. 
Ring dimensions a = l,b = 0.5. 
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Figure 2. Argand diagram for the eccentricity wavenumber with 
increasing resonant strength. Ring dimensions a = 1, 6 = 0.5. 



In between these two extremes the wavenumber behaves as 
shown in Fig. |5| which was produced by solving equation 
14UI I numerically, initially increasing before converging on 
the value of the high-strength limit. The behaviour of the 
growth rate in terms of W is shown in Fig. |3| and shows 
a maximum at a value of about W — 10, above which the 
growth rate falls and the mode stabilises at about 10^. In- 
creasing the growth rate therefore ultimately leads to a sta- 
ble mode rather than a highly unstable one. This counter- 
intuitive behaviour arises because of the cusp, which, as 
shown by Fig. deepens as the growth rate increases and 
leads to a low eccentricity at the resonant point. Since the 
eccentricity growth here is proportional to the eccentricity 
itself, a low eccentricity slows the growth of the instability 
and allows the eccentricity to be damped away in the disc. 

Whether the eccentricity is able to grow is determined 
by the parameters ab and W . Fig.|l|shows the marginal sta- 
bility curve in the plane W-ah, with values below the curve 
shown being unstable. At high W, the disc is unstable only 
for very low values of the viscosity, whereas it is unstable 
for most viscosities at lower W. 
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Figure 3. The dependence of eccentricity growth on resonance 
strength for a ring with a = l,b = 0.5, Qb = 0.1. Growth rate 
expressed in terms of the orbital frequency. 
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Figure 4. The marginal stability curve for a narrow ring with 
a = 1,, b = 0.5 in the plane W — ai^. 



3.2 Solutions of the full eccentricity 
equation-precession 

We solved the eccentricity equations 12211 . 12611 and 13311 
for their normal modes numerically using a 4th/5th order 
Runge-Kutta method with adaptive step-size. The pressure 
p and density p of the two-dimensional disc in section |5| are 
interpreted as the vertically integrated pressure P and sur- 
face density E of a three-dimensional disc. In solving the 
equations, the distributions chosen were for a steady alpha- 
disc with Kramers opacity, as described, for example, in 
Frank, King & Raine (2002) . These have the form 
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where Psc and Esc are the pressure and surface density at 
the scaling radius r^c- The tanh functions implement the 
boundary conditions that pressure and surface density are 
zero at the outer edge (they are naturally zero at the inner 
edge in this model) and contain a parameter v to allow the 



width of the drop to to be varied. In all cases, the adiabatic 
index 7 = 5/3. 

We describe the disc properties using a characteristic 
value of the ratio of the disc semi-thickness to the radius, 
H/r. The disc scale height H satisfies the relation 

P = T.n'^H^. (50) 

The ratio H/r depends weakly on radius, given by the ex- 
pression 

3/20 / 1/8 

(51) 
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where the characteristic value h depends on the system pa- 
rameters shear viscosity a, mass transfer rate M and white 
dwarf mass Afi and radius rin as 
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The constant was evaluated using a numerical solution of a 
disc with full vertical structure as described by Ogilvie and 
Dubus (2001). 

A particular binary system is characterised by its mass 
ratio q and semi-major axis a, while the disc itself is de- 
scribed by the characteristic disc semi-thickness h and the 
bulk viscosity parameter ah- The radius of the 3:1 Lindblad 
resonance is given by the expression 

'''''^(l + g)"'''^a (53) 



= 3" 



and the outer radius of the disc was evaluated using the 
prescription of Whitehurst and King (1991), in which the 
disc is truncated at 0.9 of the Roche radius as given by the 
Eggleton (1983) approximation, so 

_ 0450 

'^""'^ 0.6 + g2/3iog(l + g-i/3)- ^^^^ 

The first set of solutions obtained were for the isolated 
disc described by equation I22II . This produces a series of 
normal modes where the eccentricity propagation is affected 
only by pressure and is therefore expected to be retrograde. 
A WKB analysis of the expression was also carried out to 
act as a check on the numerical method. In fact the eigen- 
value of equation 1221 is only dependent on h and the scaled 
eccentricity distributions are the same for all values of h, 
as a transformation to a new scaled eigenvalue u —^ uj'h~^ 
completely eliminates h from the equations. 

The next solution obtained included the effect of the 
companion as described by equation 12fciL using the orbital 
parameters of the well-studied eclipsing dwarf nova OY Cari- 
nae (Wood et al 1989) which has a mass ratio q — 0.102. 
The companion is expected to produce a prograde preces- 
sion of the disc. We find that the gravitational perturbation 
dominates for the first two modes and overall precession is 
prograde, but for the remainder of the modes, the precession 
is retrograde and dominated by pressure. The reason for this 
is that, as the mode number increases, the eccentricity gra- 
dient becomes larger due to the increasing number of nodes 
and therefore the contribution of the eccentricity gradient 
term increases. This also affects the damping of eccentricity 
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Figure 5. The lowest order mode of a disc witiiout tlie resonance, 
for OY Carinae witli h = 0.003 and ttb = 0.05. 
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Figure 6. The lowest order mode of the full eccentricity equation 
for OY Carinae with /i=0.003 and ttb = 0.05. 



since this is related to the eccentricity gradient in the same 
way, and we therefore expect the lowest order mode to be 
most significant as it is the fastest growing. This lowest or- 
der mode has the form shown in Fig. |3 with the greatest 
eccentricity at the outer edge of the disc. 

Finally, we studied the full equation with terms for gen- 
eration and damping of eccentricity. In the same way as the 
previous section, the strength of the resonance can be de- 
scribed by a dimensionless parameter which we now write 
as 

W = 2gpn(rout-n„) ^^^^ 
P 

res 

For 0.05 < q < 0.3, the typical mass ratio range for SU UMa 
systems (Patterson et al (2005) propose q < 0.35) then we 
find the approximate range 20 < W < 650 with the physi- 
cally realistic value h = 0.01. The parameter W is dependent 
on h, with W oc h~^, and for the eccentricity distribution 
shown in Fig. |S] with h — 0.003 the eccentricity suppression 
at the resonant radius is correspondingly larger. The cal- 
culated range of W suggests that, in real systems, we are 
in the region where the local suppression effect of the reso- 
nance is important. The addition of the resonance to a real 
disc indeed leads to a cusp in the eccentricity distribution at 
the resonant point and the eccentricity distribution, shown 
in Fig. |S| now peaks inside the resonance at about 0.37a. 
We find the precession rate is most strongly affected by the 
mass ratio and semi-thickness, with the bulk viscosity hav- 
ing a negligible effect. 

Fig. |7| shows the effect of varying h in the disc, again 
using the OY Carinae parameters. The value of the semi- 
thickness parameter affects the relative importance of the 
gravitational and pressure terms in the expression; a higher 
semi-thickness in a hotter disc increases the effect of pres- 
sure. This dependence is linear, and for h > 0.015 the pres- 
sure contribution is large enough to make the disc precess 
retrogradely. OY Carinae has an actual precession rate of 
0.0198r2b, and therefore a low value of H/r seems to be 
more accurate, with a best-fit value h ~ 0.0035. This is sig- 
nificantly lower than the value expected from the definition 
of h in equation I5H which, for reasonable values of primary 
mass and accretion rate, is about 0.01. 



c 0.025 
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Figure 7. The dependence of precession rate on the semi- 
thickness parameter h for OY Carinae (q = 0.1). 



We then calculated the variation of precession rate with 
q for several values of h, shown in Fig. |H| with the observa- 
tional data from Patterson (2001). Once again the lower val- 
ues of h produce more accurate agreement, with h ~ 0.003 
the most consistent value, although the observational errors 
are sometimes significant. We have assumed that the effect 
of changes in the ratio rin/o are unimportant. 

Fig. |H] also shows a linear relation over most of the 
curve, with deviations at high and low q. That at high q 
is due to the effect of a resonance very close to the disc 
edge, and is strongly affected by the implementation of 
the boundary conditions in the pressure and surface den- 
sity distributions. Over most of the range, and particu- 
larly where most observed superhump systems lie, the re- 
lation is linear. A best-fit to the linear part of the curve for 
h = 0.003, the closest value to the observations, has the re- 
lation e = (0.2076 ±0.0003)g- (4.1 ±0.6) x lO-"*. This differs 
from the Patterson (2001) relation by having a small offset 
from the origin, which is due to pressure-induced retrograde 
precession of a stable eccentric mode in systems with very 
low q. An isolated system has retrograde eccentric modes, 
but in the absence of resonant excitation the viscosity damps 
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Figure 8. Precession rate as a function of mass ratio for differ- 
ent disc semi-thickness parameters. POl labels data points from 
Patterson (2001). 



these away and no long-term eccentricity develops. Other- 
wise the two relations agree reasonably well. 



3.3 Solutions of the full eccentricity 
equation-growth 

Having established the precession rates for the eccentric 
discs, we then looked at the growth rates of the eccentric 
modes in the disc. We find that this is most strongly af- 
fected by the bulk viscosity and mass ratio, although the 
semi-thickness also has some effect. Fig.|5|shows the depen- 
dence of the growth rate on the semi-thickness h for the OY 
Carinae system with a bulk viscosity Oh = 0.1. Here the 
eccentric mode is unstable at large h, but starts to become 
stable at the low values of h consistent with the precession 
rates discussed above. This behaviour is not general, how- 
ever, as shown by Fig. 1101 which shows the growth rate as 
a function of q for different mass ratios. The sharp peaks 
and troughs are due to changes in the eccentricity at the 
resonance depending on the eccentricity distribution at the 
outer edge of the disc: at certain values of the parameters 
a peak can become 'trapped' in the outer part of the disc 
and lead to a change in the eccentricity at the resonance. 
For Qb = 0.1, the general behaviour of growth rate shows 
stability over much of the parameter range for all values of 
h. 

To determine where the modes are unstable, we there- 
fore consider the dependence of the growth rate on the 
bulk viscosity, shown in Fig. 1111 This, as expected from the 
growth rate equation (I37I I. is approximately linear and gives 
an unstable mode for Qb < 0.079. To obtain a reasonable 
value of the growth, however, we need a much lower value 
of the bulk viscosity, closer to ctb — 0.01. Superhumps are 
typically observed within a day of the start of a superout- 
burst, and reach their maximum amplitude within two days 
(Kato 1996). In terms of a typical orbital period of about 2 
hours, we have a characteristic growth time {—■kQ{uj))~^ for 
an increase by a factor e which, for the highest growth rate 
in Fig. 111! (without any damping) , leads to growth by only 
a factor 1.02 over a day. This extremely low growth rate is 



Figure 9. Growth rate as a function of semi-thickness for OY 
Carinae with Qb = 0.1. 
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Figure 10. Growth rate as a function of mass ratio for different 
semi-thickness parameters. Bulk viscosity «b = 0.1. 

clearly due to the strong suppression of the eccentricity at 
the resonance. 



4 SUMMARY AND DISCUSSION 

We have shown that the effect of including the 3:1 Lindblad 
resonance in the dynamics of a two-dimensional eccentric 
disc leads to a suppression of the eccentricity distribution 
at the resonant point, which has a very important effect on 
the behaviour of the disc. The precession relation obtained, 
given by equation 1361 and replacing the dynamical equation 
can produce a much more accurate fit to the observa- 
tional data than the simple particle precession rate at the 
3:1 resonance. Numerical evaluation of the integrals in equa- 
tion 1361 1 for typical eccentricity distributions shows that the 
retrograde pressure terms are only ~ 1 per cent of the dy- 
namical term, and so the dominant effect in reducing the 
precession rate from the particle precession is the averag- 
ing over the disc rather than the retrograde precession. The 
obtained eccentricity distributions peak at about r ~ 0.37a, 
the value which Patterson (2001) suggested to match the ob- 
servational data with a dynamical precession. This appears 
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Figure 11. Growth rate as a function of bulk viscosity for OY 
Carinae with /i=0.003 



to be an explanation for the observed disc precession rates 
being so much lower than the expected dynamical precession 
at the resonance. 

However, to obtain an accurate fit to the data, we have 
to consider a value of the semi-thickness parameter h which 
is much lower than the value we would expect. It is difficult 
to reconcile h = 0.003 with the value h ~ 0.01 which we 
expect from the standard a-disc model, as this requires the 
disc to be substantially thinner and cooler. This also leads to 
the associated problem with the growth rates, since a lower 
h increases the parameter W and so decreases the growth 
rate, as shown in Fig. 1101 

It is possible that these discrepancies between expected 
parameters and those required to reproduce the observations 
are due to the simplifications involved in the model. Three- 
dimensional effects may be significant, as would be the in- 
troduction of a shear viscosity. The general behaviour of 
eccentricity will become non-linear once the linear instabil- 
ity described by Lubow (1991) has started to take hold, and 
this will modify the eccentricity distribution away from the 
linear results obtained here. In particular, non-linear effects 
will be most significant where dE/dr is large, and there- 
fore we would expect them to be important in the region of 
the cusp. This might have the effect of smoothing out the 
cusp and therefore reducing the suppression of eccentricity 
at the resonance, but further non-linear studies are required 
to investigate this. 

An estimation of the additional effects can be made 
by comparing the expressions obtained above with the non- 
linear eccentricity equation obtained by Ogilvie (2001). 
Starting from the three-dimensional, non-linear eccentric- 
ity equation obtained there and working in the linear limit, 
we use the values for the numerical constants obtained by 
numerical solution of the associated structure equations in 
that paper and obtain the following eccentricity equation 
to compare with the complex version of equation 1221 for a 
steady a disc: 



2rn 



dE _ 0.06 + 1.262i d 
'dt ~ 



dr 



Pr 



dr 



+ 



(-0.306 + 0.238i)P^ + (-0.531 + 0.938i)^ -. 

or E r 



(56) 



Equation 12211 takes the following form once the power-law 



form for the pressure has been used to calculate the deriva- 
tive; 

The effective bulk viscosity is slightly reduced from the shear 
viscosity used (0.1), a consequence of the viscous overstabil- 
ity discovered by Ogilvie (2001) which may help to interpret 
the low values of Oh required for an unstable mode. Intro- 
ducing a shear viscosity to the problem would allow the over- 
stability to take effect and might increase the growth rates 
to more realistic levels. The value of 1.262, which is the effec- 
tive adiabatic index, is shghtly reduced but not significantly 
so, which would tend to reduce the effects of pressure. Since 
we do not believe these to be very significant in determining 
the precession of the lowest order mode, it is likely that any 
effect of this change would be due to changes in the shape 
of the mode itself. 

The most significant difference between the two expres- 
sions is the term proportional to E alone, which has changed 
sign with the imaginary part of the non-linear expression 
now positive rather than negative and now describes a pro- 
grade pressure-driven precession. This term may have the 
effect of increasing the value of h which fits the observa- 
tional data, since it will tend to increase the precession rate 
of the disc. Further investigation is certainly required to in- 
vestigate this effect and determine whether a more realistic 
value of h can be obtained. There is also a completely new 
term proportional to dE/dr, which may have an entirely 
new effect on the disc dynamics. 

Other refinements to the model may lead to further 
changes in the precession rates. Numerical solutions of our 
equations with a varying outer disc radius confirm that these 
variations, which are expected to occur during superout- 
bursts, do not have a significant impact on the precession 
rate. The form of the eccentric mode in Fig.|S|suggests that 
the reason for this is the fact that the precession rate in 
equation I36I I is dominated by the part of the mode inside 
the 3:1 resonance. Variations in the surface density and pres- 
sure profiles may be more significant, and observations and 
simulations suggest that there is a build-up of mass at the 
outer edge of the disc during superoutbursts (for example 
see Horne & Marsh 1986). It is difficult to estimate the pre- 
cise effect of these variations, however, since they may effect 
the eccentricity distribution in ways which can only be de- 
termined by full solution of the eccentricity equation. 



5 CONCLUSIONS 

(i) It is possible to describe the dynamics of an eccentric 
2D accretion disc using a single equation which contains 
information about the eccentricity distribution of the disc, 
its precession and propagation, and the resonant growth and 
viscous damping of the eccentricity. 

(ii) By including the 3:1 Lindblad resonance as a delta- 
function in the eccentricity equation, we find that the res- 
onance suppresses eccentricity at the resonant point itself 
and leads to a deep cusp in the eccentricity distribution of 
the fastest-growing mode. This in turn leads to a certain 
amount of suppression of the eccentricity growth. 

(iii) By solving the eccentricity equation, it is possible to 
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obtain a more accurate fit to the observations of superliumps 
than has previously been managed using dynamical esti- 
mates alone. The closest fit, however, is obtained for values 
of the characteristic disc semi-thickness of about h — 0.003, 
much lower than the value h — 0.01 which we expect. 

(iv) Further non-linear and three-dimensional calcula- 
tions are required to obtain a more complete and accurate 
picture of the dynamics of eccentric discs in SU UMa sys- 
tems. 
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